/*
 * =====================================================================================
 *
 *       Filename:  mixbern.h
 *
 *    Description:  mixtures of multinomial distribution
 *
 *        Version:  1.0
 *        Created:  2009年06月25日 14时36分54秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Ying Wang (WY), ywang@nlpr.ia.ac.cn
 *        Company:  Institute of Automation, Chinese Academy of Sciences
 *
 * =====================================================================================
 */

#ifndef MIXBERN_H
#define MIXBERN_H

#include "ncmatrix.h"
#include "ncvector.h"
#include "kmeans.h"
class MixBern
{
public:
	MixBern(const NCmatrix<double> &mdata, int K);

public:
	void fit();
private:
	double estep();
	void mstep();
private:
	NCmatrix<double> data;
	int KK;
	int dim;
	int N;
	NCvector<double> fracs;
	NCmatrix<double> mu;
	NCmatrix<double> resp;
};

MixBern::MixBern(const NCmatrix<double> &mdata, int K)
	: data(mdata)
	, KK(K)
	, dim(data.column())
	, N(data.row())
	, fracs(KK)
	, mu(KK,dim)
	, resp(N,KK)
{
}

double MixBern::estep()
{
	int k,n,d;
	double logprob,logp,xx,mm,expsum;
	for ( n=0; n<N; n++ )
	{
		logprob = 0.;
		for ( k=0; k<KK; k++ )
		{
			expsum = 0.;
			for (logp=0., d=0; d<dim; d++ )
			{
				xx = data[n][d]; mm = mu[k][d];
				logp += (xx*log(mm) + (1-xx)*log(1-mm));	
			}
			resp[n][k] = log(fracs[k]) + logp;
			expsum += exp(resp[n][k]);
		}
		logprob += expsum;
		for ( k=0; k<KK; k++ )
		{
			resp[n][k] = exp(resp[n][k])/expsum;
		}
	}
	return logprob;
}
void MixBern::mstep()
{
	int n,k,d;
	double Nk,tmp;
	for( k=0; k<KK; k++ )
	{
		for(Nk=0., n=0; n<N; n++ )
		{
			Nk += resp[n][k];
		}
		fracs[k] = Nk/N;
		for( d=0; d<dim; d++ )
		{
			for( tmp=0., n=0; n<N; n++ )
			{
				tmp += resp[n][k] * data[n][d];
			}
			mu[k][d] = tmp/Nk;
		}
	}
}
void MixBern::fit()
{
	for( int k=0; k<KK; k++ )
	{
		fracs[k] = 1./KK;
	}
	kmeans(data,KK, mu);

	
}
#endif
